NASA TECHNICAL 


632 NASA TT F-670 
TRANSLATION ur. 


NASA TT F-670 


PROBLEM SOLUTION BY 
THE “LARGE-PARTICLE” METHOD 
by K. A. Vedyashkina, Z. F. Levina, S. P. Lomnev, 


G. P. Prudkovskiy, 1. V. Rastopchina, G. V. Ruben, 
and V. V. Yurchenko 


Transactions of tbe Computer Center 
Academy of Sciences, USSR 
Moscow, 1970 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION * WASHINGTON, -D. C. € JULY 1971 


PROBLEM SOLUTION BY THE 


"LARGE-PARTICLE" METHOD 


By K. A. Vedyashkina, Z. F. Levina, S. P. Lomnev, 
G. P. Prudkovskiy, T. V. Rastopchina, 
G. V. Ruben, and V. V. Yurchenko 


Translation of "Resheniye Zadach Metodom Krupnykh Chastits " 
Transactions of the Computer Center, Academy of Sciences, USSR 
Moscow, 1970 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the National Technical Information Service, Springfield, Virginia 22151 
$3.00 


PROBLEM SOLUTION BY THE 
"LARGE-PARTICLE" METHOD 


Editor-in-chief: S. P. Lomnev, 


Doctor of Physics and 
Mathematical Sciences 


TABLE OF CONTENTS 


intToductiobi ias iaa RC ak aa 117177” 1 
The Movement of Charged Particles in an Electric Field ............ 3 
Two-dimensional Problem of the Behavior of a Cloud of 

Charged Particles ina Magnetron ........................ 6 
Investigation of the Particle-in-Cell Method .................... 17 


Some Questions in the Calculation of the Evolution of Stars .......... 51 


INTRODUCTION — 


Methods, based on replacing the continuity equations or their analogs by "large- 
particle" equations of motion, are now finding application in the solution of problems in 
electrodynamics, gas dynamics, hydrodynamics and even solid-state physics. This 
method has given good results in a number of problems where classical methods en- 
counter difficulties. 


Two steps can be identified: 

1) the "particle" motion during At; is calculated in a given field; 

2) from the new "particle" configuration the fields are calculated for a new deter- 
mination of the motion during At, +1 


The "large particle" itself has various geometrical shapes: point, ring, disk, rod, 
filament. 


The details of the problems and the necessity of adapting different variations of the 
method to them give birth to an abundance of procedures and a confusion in terminology 
is possible. It seems appropriate to us to identify three basic methods in terms of the 
nature of the interaction to be calculated: 

1. The straightforward "large-particle" method. The fields are known functions of 
the coordinates and velocities. The problem is solved in one step: only the particle 
motion is calculated, Because of its simplicity this method has been used to solve sev- 
eral electrodynamic problems. 

2. The "average density" method. A region is marked off into a three-dimensional 
lattice; the density is calculated for an arrangement of a finite number of test particles, 
whose motion in the time interval At; E is found from the dynamic equations with the 


forces determined by this density at a previous instant of time b. 


3. In the Harlow "particles in cells", or PIC, method the solution is carried out in 
three steps: 
a) preliminary values of the fields are found during the time interval At, in an /4 


Eulerian network of cells, fixed with respect to the observer; 

b) the motion of the masses is determined from the equations of motion in a 
Lagrangian network of particles; 

c) the field corrections, associated with particle motion, are calculated. 


* Numbers in the margin indicate pagination in the foreign text. 


The PIC method has been found to be quite consistent in a number of gas-dynamics 
problems. 


This collection of papers contains examples of problems solved by these three 
methods. They do not represent a detailed physical investigation, but, instead, they 
illustrate the calculational possibilities. 


THE MOVEMENT OF CHARGED PARTICLES /5 
IN AN ELECTRIC FIELD 


K. A. Vedyashkina, T. V. Rastopchina 


ABSTRACT: Application of the PIC method to calculation 
of motion of charged particles in an electric field is ex- 
amined. 


In some problems it is advantageous to use the straightforward "large-particle" 
method. The volume V, containing charges at t = 0 is divided into small volume elements 
AV (the volume elements can be equal in size), and each charge, contained in AV mn 


is concentrated at its center, It is considered as a new "enlarged particle". Altogether, 
M "particles" are obtained, whose motion is determined by the systems of motion 
equations: 


dP > R.H 
de ^ ‘mjEn* ... 
E m Eon +E Im? (1) 
H,, s las H ym? m = 1. 2, , M. 
Eom? Hom are the external electric and magnetic fields at the Ry th point; Elm? Him 


are terms that take account of the interaction of the m-th "particle" with all the others. 


Averages over the particles within AV (m — 1, 2, ..., M) are taken as the initial 


conditions. In this it is assumed that during the time 0 « t « T the active particles in the 
volumes AV ny act as a unit; this is a permissible assumption for finite T (the particle 


transit time) and for sufficiently small AV when long-range interactions predominate 
over short-range within AV ra The permissible size of AV depends on the coefficients /6 
and on the spread of the initial velocities in AV. mə It is determined by the convergence 


of the solutions of Eqs. (1) for an additional decrease in AV [1]. 


The effect of near collisions of the i-th and j-th particles is excluded by the 
introduction of a "protective cube", the dimensions of which are determined on the basis 
of an experimental analysis or the following simple considerations. Let the scattering 
of a particle with charge x. by particles with charge e: be expressed by the well-known 


Rutherford formula and the scattering force is 


r 
hə [f e yı: e itj 
2.297 T 2° 
D q +? m, 


for r, «<r we obtain 
p2,p? R r? "1 max’ .... 
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r?. 
i 
a "protective cube" size of the order of Ti [2]. 


As an example let us consider a constant-current electron beam, moving in space 
with a uniform electric field Eo = const, Eo =0, Eso = 0 for zs L. It is simu- 


lated by "large particles", emerging from the slit in time At, whose charge is 
4 =] (t,) At, Vio (tJ = Vo; t “to + (k-1) At, Z 10 (t =Zo- 


Let us assume that at the pu instant of time one "particle! emerges, When Zn.” L, 


the "particle" no longer needs to be considered, Then the equation of motion will be 


: i 

Zm= AglE + Ej, E: im” : iz ni (2) 
i fm mi 
O A kel 2, see K, 


where A o -—— : e is the electron charge and m is the electron mass; c is the velocity 
c mg 
of light; 0 = z(t) = L; OL is the field-filled gap; differentiation is with respect to ct. 
In the process of solving (2) the particles are accumulated in the interval 0L; At 


is selected from trial calculations, Each particle is surrounded by an &-region to 
exclude collisions. 


Graph of p (z): 


l-forI-1,,2-forl 61 


0? 0? 3 - for I= 91). 


To calculate the charge density p with 1% accuracy it was found adequate that 
M = 300, Act “8, vg was taken equal to 0,01c and L = 10 cm, The figure shows the 


distribution p(z) at some instant of time. It is not hard to see that as the current in- 
creases, p(z) becomes more and more nonuniform and fluctuations appear, caused by 
Coulombic interaction. 
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TWO-DIMENSIONAL PROBLEM OF THE BEHAVIOR OF A 
CLOUD OF CHARGED PARTICLES IN A MAGNETRON 


K. A. Vedyashkina, G P. Prudkovskiy 
ABSTRACT: Application of the PIC method to behavior 
of electron beams is discussed. 


Magnetron design by computer is now possible (see [1]. This paper considers the 
self-consistent motion of a large number (of the order of a thousand) of charged particles, 
representing electron beams. 


In this paper the analysis is applied to the nigotron type of continuous-operation 
oscillator [2, 3]. 


The interaction region is a rectangle which is subdivided by a constant-interval 
gridwork for calculating the space-charge field by the "average-density" method [4]. 


The charge density p (x, y, t) is calculated from the disposition of test particles at 
the end of a time interval At and determines the field acting on the particles during the 
subsequent interval. 


The equations of motion can be written in the following form: 


y -AQCEL- XH, (1) 


the time is in length units (ct instead of t). 


Let us examine E = Es + E A” E p* We assume the d-c field is uniform Es = 


u 
Foy =~ EE , where ua is the anode potential and d is the width of the active space. 


Such an approximation is permissible for a one-row negotron in which the cathode is 
continuous, as in the normal magnetrons. 
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In the high-frequency field let us isolate one spatial harmonic (see [2] ) 


F,.- E, shaveosgtx —Uppt): | (2) 


E, E, chay sinqix Uf). 


where g = = ,V.,7c/A, f=1 cm is the period of the system. 


ph 
The space-charge field E p is determined by the potential, for which 
AD(X, v, t) - -Anp(x,y, 0), (3) 
with 9 = O for y “0 and y =d, 
(x, yət) = 9 (x+1l, y, ty (4) 
The magnetic field is uniform (H, — H). 


The portion of the cathode emitting into the nigotron is small: ó xb (bis the 
system dimension along the z-axis). In each time interval k particles, equidistant from 
the x-axis over the segment 6 (6=0.2727, b — 107), are emitted simultaneously. 

Let us examine possible ways of solving the equations of motion, The method 
recommended in [1] uses à special form of Eqs. (1). The Eqs. (1) are written in 
complex form 


z«iw,z -AgE, (5) 


where MW. de A¿H; Z=X>+iy; E = E. + E; 


we find 
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Assuming E(£ ) = E(t p! (p is the interval number) in At, we obtain 


y e 
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Besides the analysis in accordance with the scheme given, the equations were also 
solved by the Runge-Kutta method. 


Let us compare the errors in one interval At of the analysis. According to the 
method of [1] 


A -O9,.w. At 
o dE c t? 
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Thus, both errors are of the order of At”, Sample calculations showed that both 
procedures require the same amount of time to calculate trajectories with the same 
interval, The interval, required for the calculation, was selected from a comparison 
of the results; it must not be greater than T,/50. For the calculations it is recommended 


that an interval be taken such that a particle passed through, on the average, one cell of 
the point lattice. 


The calculation of the space-charge field is based on solving Poisson's equation in 
finite differences. For the i, j lattice element we have 


2515-55 ae ae me AC j-1 =—{-- 9 
22 p-n... (9) 
where 7 = d > ho zs. pi = i. : M is the number of elements along the 
N 3 M . > j 4hTb , g 


x-axis; N is the number of El eure along the y-axis; m. . is the number of particles 


in four rectangles surrounding the point (i, j). The 2. of the test particle is chosen 
from the condition 


l,^t = qk. (10) 
The boundary conditions are 
Tope COLLO I n (11) 


The field is determined by solving a set of N x M linear equations. In [1] the 
potential is represented in Fourier series. In our work we used the Seidel iteration 
method, which is suitable for this problem because the values öl E. obtained in the 


previous interval and differing from the desired values, are ‘akon as the zero approxi- 
mation, Therefore the iteration process converges quickly. The field components can 
be found from the formulas 


ACIE E (12) 


—-4$. 
E "ede Tj h^). TERT hz 
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The calculations were performed on a BESM-6 computer. To ascertain the param- 
eters required for the analysis the trajectories were calculated for different numbers of 
particles and different numbers of lattice elements. It is found that a 15x 30 lattice is 
adequate for a current I, = I, and u, = 20 KV; at least 6 particles must be emitted during 


the interval At = T 760. For a higher current, as seen from Fig. 1, thereisan 
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appreciable difference between the trajectories plotted for 15x 30 and 30x 60 lattices. /16 
At 5 amperes the character of the trajectory is altered qualitatively when one goes over 
to the more precise solution, 


It should be noted that the effect of the space charge is a slowing of the drift of 
particles to the anode. Characteristically its effect is to reduce the amplitude of the 
cyclotron oscillations, In Fig, 2 the particle coordinates are shown for a fixed instant 
of time. They form the typical magnetron tongue, Figure 3 shows the tongue for a cur- 
rent that is four times higher than in Fig. 2. The tongue "blowing-away" effect, described 
in [2], is clearly seen, Figures 4, 5 and 6 illustrate the change in anode and cathode 
current for different analysis parameters and for different emission currents. Figure 4 
shows current graphs ignoring the space charge for 6 and 12 particles emitted during At, 
An increase in the number of particles causes a flattening of the graphs. It is significant, 
however, that a simple averaging of the oscillations gives the same effect. 


In Figs. 5 and 6 the currents are 3 A and 5 A, respectively; graphs are given for 
different analysis parameters. As the analysis is refined, a splitting of the cathode- 
current maximum is observed, The calculation time amounts to 45 minutes (30 minutes 
when interaction is ignored) for 10 T Periods (600 time intervals) for an emission of 


12 particles/At and a 30x 60 lattice. 


45 45 


A JNE Y . 
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45 55 65 


Figure 5. The Dependence of LA, on t for pc 3 A for 


k — 6, 15x30 lattice; k € 6, 30x 60 lattice; k = 12, 30x 60 
lattice. 


This paper is methodological in nature. The calculations in accordance with the 
program formulated indicated the possibility of studying the physical processes, 
occurring in a magnetron-type oscillator, with the aid of a BESM-6 type computer. 
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45 | 55 6,5 15 tí, 


45 55 — 65. 75 tj, 


Figure 6. The Dependence of L/L, on t for 1. =5A for 


k=6, 15x30 lattice; k = 6, 30x60 lattice; k = 12, 30x60 
lattice; k = 24, 30x60 lattice. 


The equations of motion were integrated by the Runge-Kutta method, the Poisson equa- 
tions — by the net-point method using the Seidel iteration method. The necessary analysis 
parameters are calculated. Satisfactory results are obtained for the simultaneous motion 
of 2000 particles and a calculation of the field from a 30x60 lattice of points. Doubling /18 
the number of particles makes it possible to explain the more subtle effects in the current 
relationships. Calculations show that a consideration of the higher harmonics of the d-c 
and high-frequency field does not increase the volume of the analysis too drastically. 

This makes it possible to proceed to a study of actual double-row nigotron designs. 
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INVESTIGATION OF THE PARTICLE-IN-CELL METHOD 


V. V. Yurchenko 


ABSTRACT: This paper discusses light computational schemes 
of the PIC method, as applied to supersonic flow of a gas. 


1, INTRODUCTION 


In recent years, many papers have been written devoted to many variants of the 
method of "large particles" [1-5], in particular, for the solution of hydrodynamic prob- 
lems by the particle-in-cell method [1, 4-8]. This method was tested by F. H. Harlow 
in 1957 [1], was published in a final edition in 1963 [7], and was further developed by 
co-workers at the Los Alamos Scientific Laboratory. During these years, he suc- 
ceeded in solving numerous interesting physical problems. However, a detailed des- 
cription of the method has not been openly published. Therefore, attempts to repeat 
the solution have encountered serious difficulties. The necessity has arisen for a 
more painstaking investigation of the particle-in-cell method in order to clarify its 
theoretical details. 


2, DESCRIPTION OF THE METHOD 


The particle-in-cell method is connected in a natural manner with various approaches 


to the solution of the hydrodynamic problems of Lagrange and Euler. 


As we know, in Lagrange's representation, a liquid is partitioned at the initial in- 
stant of time into finite zones. Deformation of the network of zones in the course of time 
is connected with the motion of the liquid, which depends on a given numerical method 
for solving the differential equations of hydrodynamics. This technique exhibits a high 
degree of accuracy and is convenient for the solution of nonstationary problems. It 
enables one to follow the motion up to and including the fine structures of the flow. How- 
ever, this approach has its defects, the most significant of which is its inapplicability in 
cases in which the medium undergoes strong deformations, large relative displacements, 
or is subject to other factors leading to a strong distortion of the original network. 


In Euler's representation, the liquid flows through a network of finite-difference 
cells that is fixed relative to a stationary (laboratory) system of coordinates. An im- 
portant advantage of Euler's approach is the ability to calculate the motion of the liquid 
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under strong deformations and large relative displacements. But even this technique has 
a number of defects such as the difficulty of determining the fine structure of the flow, 
the difficulty of calculating the interfaces of the media, etc. 


The particle-in-cell method is one of the attempts at combining these two different 
approaches, keeping the best characteristics of each, In calculations made with this 
method, the region of the solution is partitioned into Eulerian cells that are stationary 
with respect to the observer. One introduces a Lagrangian network of particles repre- 
senting the elements of the fluid moving through the Eulerian network, The Eulerian 
network then describes the variable fields of flow and the Lagrangian network describes 
the parameters of the fluid itself, Obviously, the use of the two networks presents high 
requirements on the operating speed and memory of a computer. One should note the 
great graphic obviousness of the method, whichis close to an analog y between the calcula- 
tion and actual experiment (we impose initial and boundary conditions and then need only 
observe the development of the process). After each cycle, the result determines the 
real flow of the fluid at the corresponding instant of time. This fact points out the broad 
possibilities of the particle-in-cell method, especially in the solution of nonstationary 
problems with strong deformations of the flow. 


In calculations made by the particle-in-cell method, the medium is represented in 
the form of discrete point masses, that is, particles moving through the Eulerian network, 
Here, the mass of an individual particle does not change. This makes it possible to elimi- 
nate from the complete system of hydrodynamic differential equations the equation of con- 
tinuity, expressing the law of conservation of mass. Calculation of the change in the 
pattern of the motion of the medium involves time steps, or cycles. A cycle occupies a 
time interval ót. Each cycle is broken into two stages. The first stage involves the cal- 
culation of the flow of the fluid in Euler's representation. In this stage, the position of 
all particles is considered completely fixed. Therefore, in a system of equations de- /21_ 
scribing the motion of an ideal fluid, the terms reflecting the displacement of the medium 
are omitted, and the system takes the following form 


gradP; .. ə Ediv V, (1) 


where I is the specific internal energy of the medium, The result of the calculations at 
the Eulerian stage is the determination of the preliminary new values of Y and I at the 
end of a cycle. At the Lagrangian stage, one considers the motion of particles during 
the time Öt of the same cycle. This stage takes into account the displace ment of the 
medium, In the course of the motion, the particles transport their energy and momentum 
from one point in space to another. As a result of this redistribution of the masses, 
energies, and momenta, there will be a new change in the fields of the velocity and 
internal energy. The new values obtained for the fields of flow and the parameters of 

the fluid will describe a motion of the fluid at an instant of time corresponding to the 

end of the cycle. 


At first glance, it may appear that the particle-in-cell method does not yield a stable 
solution since any solution of the system (1) other than that of undisturbed flow is ex- 
plicitly unstable. However, as was shown in [8], the redistribution of the energy and 
momentum at the Lagrangian stage is equivalent to the addition of a number of terms in 
the equations of the system (1), some of these terms corresponding to the effective 
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viscosity and thermal conductivity, This explains the stability of the solutions obtained 


by the method. 


3. THE BASIC COMPUTATIONAL FORMULAS 


In the present article, the particle-in-cell method is applied to a specific problem 
of supersonic flow of an ideal gas around a cylinder, the axis of symmetry of which is 


directed along the incident flow, The region of the solution (see Fig. 1) is partitioned by 


means of a rectangular Eulerian network into IxJ cells, Since the problem is axi- 
symmetric, it reduces to the two-dimensional case, and the cells of the Eulerian net- 
work are annuli of rectangular cross-section. Each cel (i, j) is characterized by the 


following variables: 


u. .andv. . 
1, ] 1, ] 
M. 
1, ) 
L. 
1, ] 
P. 
1, ] 
E. 
1, ] 
Kz. . and Kr. 
> J 1, 


For each cycle, the initial conditions are the distribution of the particles and the 


values of the quantities M. 


1 


the components of the rate of flow in the directions z and r; 


the mass of the gas 


the specific internal energy of the gas; 


the pressure; 


in a cell; 


the total energy of the gas; 


the components of the momentum in the directions z and r. 
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Figure 1. The Region of the Solu- 


tion. 
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preceding cycle or are given as initial conditions 


.. Which are either determined in the 


of the problem, 
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The system of equations (1) that describes the Eulerian stage of the calculation 
becomes, in cylindrical coordinates 
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where P. .is determined from the equation of state and Vi i is the volume of the 


(i, j) cell. 
When one specifies the boundary conditions for P, u, and v, one completely deter- 
mines the system (3) - (5). The solution of this system yields u; ; v ; and I, i 


These are the tentative new values of the corresponding quantities, The boundary condi- 
tions are determined by the laws of conservation of the total energy of the system and 
the momentum. The question of the choice of boundary conditions will be examined in 
detail below. The Eulerian stage in the calculation ends with this. 


Before turning to the Lagrangian stage, we need to carry out some supplementary 
operations. Specifically, we need to calculate the tentative total energy and momentum 


for each cell 


-2 -2 

Uu . e +> V. . 
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At the Lagrangian stage, the displacement of particles is calculated in accordance 
with the values of u, , and v, found during the cycle time ót. As a result of this 
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123. 


displacement, the coordinates of the particles change. Here, a particle can either re- 
main in the same cell from which it began its motion, or shift to another cell, or com- 
pletely leave the region of solution through the open boundary. Each k-particle that begins 
its motion in the cell (i, j) takes with it a mass m, and (in proportion to that mass), an 


amount of momentum and energy: 


Depending on where the particle is as a result of the displacement, the quantities Kz, op, 
Er, ; Es, . and M, j do not change if the particles do not cross the boundaries of the 

cell (i, j) and they 1 by an amount (7) if they do (in which case M; j changes /24 
accordingly). If the particle is in a new cell, then the quantities (7) and m, are added to 


the corresponding quantities of that cell; if the particle leaves the region of solution, 
then (7) and m, are "removed" from the system. With this redistribution of mass, 


energy, and momentum, the total change in them for the entire system is exactly equal 
to their flow through the open boundaries of the region of the solution. Integration of the 
motion of the particles can be done in various ways, which will be examined below. 


As a result of the displacement of the particles of mass, the total energy and momen- 
tum in the cells assume, in general, new values 


M^. E’., Kz.. Kr. .. 


5)? 55] 


These values completely determine the field of flow at the end of a given cycle: 


: Kz; ; sti Kr, j | 
8:1) Mo by] M; (8) 
M E _ ud” : o 
ij Mi 
Thus, the values Mi E. "ui ; vi, j and IH E obtained as a result of calculations in 


the Lagrangian stage E MER the solution of the problem for the corresponding instant 
of time and they serve as initial conditions for the following cycle. 


4. THE CHOICE OF BOUNDARY CONDITIONS 


To specify the boundary conditions for the system (3) - (5), it is necessary to surround 
the entire boundary of the region of the solution with a number of fictitious cells as indi- 
cated in Fig. 1. The specification of P, u, and v in them will then completely determine 
the system. 
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Our region of solution has, as one can see from Fig. 1, the following boundaries, 
each with a differing in physical meaning: 


OR, the incident flow, represented by the fictitious cells (0, j), where i = 1, 
e... 3 J; 


RB and Bj... the open boundary, represented respectively by the cells (i, J + 1), 


where i=1, ..., I, and (14 1, j), where j=j, * 1, ES 5 


Oi,, the axis of symmetry, represented by the cells (i, 0), wherei-1, ..., iyə 


Ena 
Ai, the vertical boundary of the body, represented by the cells Go + 1, j), where 
)9 155 UE 
Alp the horizontal boundary of the body, represented by the cells (i, İr), where 


i=ip+1, ... y L 


Furthermore, during stabilization of the gas flow, cells may be formed that are not 
filled with the medium. The determination of the quantities P, u, and v for such empty 
cells also calls for preliminary investigations. 


AS one can see from [3, 7, 8], we can specify conditions on the boundary and in an 
empty cell in various ways. One combination or another was used in the above analyses 
of the physical phenomena occurring in the flow at the boundaries of the region and in the 
vicinity of an empty cell. Then, the momentum and energy of the entire system were 
investigated, The results are given in [7, 8]. However, in these last communications, 
there remains some arbitrariness and contradiction in the choice of boundary conditions. 


As a first approximation, we can obtain the conditions from an analysis of the 
momentum of the system and the total energy flow in the gas through the boundary of a 


cell, The change in the component of the momentum along the axis and the total energy 
of the gas in the cell (i, j) during time ót are then as follows: 


When we use equations (3) - (5), these equations take the forms 


- V, “St T TE 9 
5Kz, 7 Vj gaz Ciao "gen" ic 
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uc aos 
6t 1 by] ta] rı : 
di Vie; AZ T 2 y Pj su; ) T (10) 
Dos +t o Y 
6t ] (5 5] "a AS 
t Vij Ar pu 2707 Po, "WA. 231- 4 
U. . -U; . 
; St 1 Ed MODAL DRA De ps, TERATE 
B EY İB 5” Ps ut 21— 14] 


The expression (10) can be written 
^E; ; < Ey Pye Fin Frye 


where Es — E 


rv 2re the total energy flows through the boundaries of the cell (i, j) 
(see Fig. 2): 


Fy 27: le x NEN "T (11a) 
tus = 5 Yd UZ ^] 5—- t E 29 (11b) 
F 1 Vj ona mu ( - qa ıl (11c) 
Fry = - O jay eg eh uni 2) : (11d) 


Figure 2. The Cell 
(i, j) in the Eulerian 
Network, 
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Since the boundary values of P for equations (3) and (4) and equation (5) can in 
general be different, we need to note that the pressure in Eq. (9) and in the left-hand 
terms of the expressions (11a) - (11d) is taken from the equation of motion (3) and that 
in the right-hand terms of (11a) - (11d) is taken from the energy equation (5). 


Let us now look at the jth series of k cells (i, j),..., (İL: j) and let us assume that 
cell (i, j) of this series is empty. Then, the total change in the momentum in these cells 
during time öt will be given by 


5t b Uii 
s ETA Y POLEA WES o o! n 
AKz, = aK z; ; V (12) 


B=, 


p. SP. 
gti "ipio Dea 7 Trag 
) 2 m 


do 


The first two terms in (12) are equal to the impulse imparted to the medium in the above 
series of cells by the forces applied to the boundaries of that series, The last term is 
the undesirable extra impulse that the empty cell imparts to the medium. It vanishes 


under the condition 


_p 
ae 5” ° (13) 


Let us now look at two possible variations of the boundary conditions for different bound- 
aries and an empty cell. 


1. An empty cell. Suppose that cell (i, j) is empty. Since M, j = 0, the equations 


corresponding to that empty cell vanish from the system (3) - (5). The effect of the empty 


cell on the system can be eliminated by satisfying condition (13) and assuming that the 
energy flows into the empty cell from neighboring cells are equal to zero. This last 
condition is satisfied when the flows (11a) - (11d) for the cells adjacent to the empty one 
are equal to zero. Thus, we obtain the following system: 


Pup Ej á 
= Uu; l,j “7 0 
by] 2 T -b/ ij : 

u. - + Hj j 

i+1,/ t+1,/ : 
2 əə ep ag 

7 V 
Poy iş ğal SIE SP; jet "hj ] x 2 - 0 

: i i 2] + 1 

P; : + UY 

; J—1 (,/—l 2 " 
e t PL ey mC + = =f), 
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par. 


The solutions of this system yield all possible conditions eliminating the effect of the 
empty cell under the requirements mentioned above, These necessary conditions are 
as follows: for a solution of equation (3), it is necessary that 


Po. -0, P - A, (14) 


i-1,j A2 


where A is an arbitrary constant that can, in particular, be zero. Equation (4) can be 
solved under the unique condition 


P. >» =): (15) 


To solve equation (5), we need to satisfy only one of the following four combinations of 
conditions: 


Po, TT P, s; 2p nion 255755: (16a, b) 


1) 
5—. 7. — 


TT, 


2. The incident flow, The boundary conditions for the fictitious cell (0, j) can be 
obtained from two considerations, namely, the constancy of the parameters of the /28 
incident flow or the constancy of the momentum imparted by that cell to the system and 
the constancy of energy flow through the boundary. In the first case, the boundary con- 
ditions for equations (3) and (5) have, respectively, the following forms: 


Po EP: Mo j=Voo. (17) 


In the second case, they can be obtained from the conditions of constancy of the momentum 
from (12) and the energy flow (11a): 


Uu, r-u R 
t l.] 1,/ 5t 
2j ao, ete dü 1o) "Vj az Pe Vut 


Solution of this system yields the boundary conditions for equations (3) and (5), respec- 
tively: 


IS E (182) 
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u ; + M 1 P u - +4 of 
-—- 1, of so a ze l,j Le} ) (18b) 


3. The open boundary. Let us look at the case i =1, m analogy with vhat was 
done above, we can obtain the boundary conditions for equations (3) and (5) in the 
fictitious cell (1 + 1, 3) from an analysis of (12) and (11b). When we equate the second 
term in (12) to the momentum imparted by the external forces and equate the flov of 
energy (11b) through the boundary i = I to the flow of energy in the cell (I+ 1, j), we ob- 
tain the folloving system of equations for determining the boundary values of P and u: 


y st 50 "T 
Jaz 2 I+1) $ 


u - +H 
1 5t lf lyf 
DT 2 «B, ja Vii AZ Pron; gu 


From this we obtain the desired boundary conditions for equations (3) and (5), respectively: 


Pu “Pi? (19a) 
1 HE Bhi (19b) 


For the open boundary i = J, the flow of energy (11d) through that boundary can, in 
analogy with what was done above, be equated to the flow of energy in the cell (i, J +1), 
which yields the equation 


b. 2 V. \ \ 
- i, f i, f A i 
eb. 1 at), jy v 01 (8 SES). 


= - AS p . Ar 
= 2sAIUE yiv, (RSF). 


The solution of this equation yields the boundary conditions for equation (5) on the 


boundary i = J: 
-1 
VD. , +1, 
fn: 2 dM 
2 YN ZU (20) 
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where P J41 


are not determined, 


are the boundary conditions for equation (4); in the present case, they 


4, Boundaries of the body. Let us look first at the vertical boundary i = Ine In 


this case, the boundary conditions are taken under the assumption that the flow does 
impart or remove heat from the body. This means that the energy flow (llb) is equal to 
zero. Thus, we have the following equations 


P. . 4 P. : 
T 5t bp] iy tle) | | . 
J AZ 2 LAZ sI? 


Hr T 7 inə / du. 271 


From this we get the boundary conditions for equations (3) and (5), respectively: 


Up etj D Pa. sj (21a) 


TUS type 22b 
MES E : : i 


To determine the conditions on the horizontal boundary of the body j = İr from the same 


considerations, we equate the flow (llc) for the cell (i, İr + 1) to zero and obtain 


U: - RU: 
P.. qu UA IN T | 2 .0 
ty} I 2 DL .] TESI " Zin " EA 


From this equation, the boundary conditions for the equation (5) on a horizontal rigid 
boundary assume the form 


Ue. 77 
lsjart I ba +1! bey 
OT ibi dt PONE E PL EN (23) 
ı N “) r B Y C, | r 
where the P, , are the boundary conditions for equation (4). Their determination re- 


i, j 
Aa 
quires the use of supplementary considerations. 
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5. The axis of symmetry, On the boundary j= 0, the flow behaves just as in the 
case of a rigid horizontal wall, Therefore, we can very simply use the results of section 
4, which yield as boundary conditions for equation (5) on the axis of symmetry 


: KU. 
Vo 5 (24) 


where the P. are the boundary conditions for equation (4), which require additional 


1, 0 
investigation, 


Thus, we have found all the necessary boundary conditions for the system of equa- 
tions (3) - (5). It should be noted that in the approximation of the equations of hydro- 
dynamics by means of our system of finite-difference equations (3) - (5) vvith the boundary 
conditions (14), (18a), (19a), and (21a), the change in the momentum of the medium of the 
region in question is exactly equal to the momentum imparted by the forces acting on the 
boundaries of the region (we are considering the momentum and the change in momentum 
along the z-axis), This total change in momentum can be obtained by summing (12) vith 
respect to j, taking into account the external forces that we specified a priori, when 
deriving the boundary conditions. This change will have the form 


/ / Et 
AK, «P, s]? (an? kəf: rk Hin. (25) 


As we can see, the expression depends only on the external forces. 


It follows from the law of conservation of energy that the change in the energy of a 
fluid in the region of solution under consideration during this time must be equal to the 
flow of energy through the boundary of the region. It is necessary to check whether or 
not this condition is satisfied for our approximating system of equations (3) - (5), with 
our boundary conditions. To investigate this, let us look at the three nonempty cells 
(i, ), (1-1, j), and (i, j +1). Obviously, this flow of energy from the cell (i - 1, 1) into 
the cell (i, j) must be equal to the flow of energy into the cell (i, j) from the cell (i - 1, j). 
This holds also for the cells (i, j), and (1,j +1). These two conditions can be written in 
the following form: 


F, ‘yj x PH, id] : Eu ini ¡+ a iy / ş (26) 
where F4 i3 is the flow Fy (see Fig. 2) for the cell (i, j). If conditions (26) are satisfied, 
then, in the summation of ‘AE, . (10) over the whole region (which does not contain /31 


empty cells), the "— flows will cancel each other pair-wise and there will 
remain only flows of energy into the region through the external boundary; that is, the 
law of conservation of total energy will be satisfied, Conditions (26) for the system of 
equations (3) - (5) have the following forms: 
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NN (272) 


Ü. . NULL 
Pig Y Va” E 5i VSA T7 Y/7 


ve (4) 


p? "ni lif, AN pl far 
“bi 2 4j zi dini” E 


(27b) 


where po, , for example, denotes the pressure ascribed to cell (i, j) when equation (3) 


is solved. This partition is necessary since it is possible to ascribe various values of 
the parameters P, u, and v to a cell during the solution of the individual equations of the 
system (3) - (5) (as was done in the derivation of conditions (14) and (16a) - (16d)). 


Equations (27a) and (27b) are valid, as one can easily see if and only if 


(3) (5) (3) ( 4) (5) (4) (5) 
yi P. du. Q 


i—1,] m i-1,J : ts] 191 ij > işğal T Jl 
- cue : ün. 
( 5) | by} iy} : ( 3) t—1,] i-1,/ : 


”, 
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—. 
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” 


Since (27a) and (27b) are written for any three cells, the preceding equations can be 
generalized as follows: 


(3) (4) f 5) 
X = ae, P; , (28a) 
& ES o (4) (28b) 
od 2. 25 “0: ı (5: Le Vo + Vi) 
CE CS t.) EE au 


Conditions (28a) and (28b) mean that, for exact conservation of the total energy, the — /32 
parameters P, u, and v in an arbitrary interior cell of the region must correspond 
uniquely to that cell, Therefore, henceforth we shall not distinguish between them, 1t 
follows from condition (28b) that, in approximating the energy equation by means of the 
finite-difference equation (5), the law of conservation of energy cannot be satisfied. 
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Equation (5) must be replaced with 


es 


bə - lij bj p, Mir 7 Mia.) 
5t i) 5] 242 
^ : 1 ^ F 3 (29) 
7.0 MES MA Loe Ula m 
2 l : 
2ar lj - z) 
where 
a "m $ “7 . ^ U jj * Vij (30) 
ute en eh uu t co 
by J 2 5] 2 


, tt ^ ^ 
Fig 73 a "ü,i in Pág Mah 912) 
b ^ m 
Fy ij =~ 3 Vi as I cA 2. Hd” HELL (31b) 
Le A cd ^ ; 3| 

Fn ij 72 -ATAZ CLE 0; (i3) 7 VU, j-1 (j- a). (310) 

F 1 ATAZSİİD U pt +P. (V . jrl . 31 
iya” ome ILE LAT 9] m pul”? (31d) 


We now find the boundary conditions for equations (3), (4), and (29) under which the 
momentum and total energy of the system are conserved, 


6. An empty cell, Suppose that cell (i, j) is empty. Then, after summation of 
AE, j over the entire region, the energy flows between the interior cells will cancel 


each other pairwise. In the empty cells (i, j) the flows of energy into the empty cell, 
which in the summation yield 


Fry i-)j ' Fi m 1/ T Marigot 1 Piv i,i-i ° 


remain uncompensated, If we make this fictitious energy flow equal to zero, we obtain 
the following equation for the conditions in an empty cell: 
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/ dane 0: T D > 
ez rə 2757 Honi ) m siy t TEM 177 yo” 
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When we take account of Eq. (13) in Eq. (32), we obtain the following two variants of 
conditions in an empty cell when the total energy and momentum are conserved: 


P;j - 0; Pia E P; aj , Pj. po” Bi (632) 
Here, a, and A ; can be arbitrary; 
P 5j 0; ee ~ ZEE "ad 2 (33b) 


^ . 
Here, u, ¿“an be arbitrary. 


If, in addition to condition (32), we require that the flow of energy into the empty 
cell from each adjacent cell be equal to zero, this imposes the supplementary require- 
ments on the parameters of the flow in the empty cell: 


Piu * 
ij Eid NO iy - 0; 
^ ^ 
P, iiij T i+1,/ "ig =): 
A A 
vw. : 2 
5 Vi j-1 + Puja 0 7 1 + 2j E z) = 


PET v 2 
ls] UL. Mur Ui j ( ES - Ü, 


(33c) 


(33d) 


where A is an arbitrary number, 
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7, The incident flow, The boundary conditions in a fictitious cell (0, j) are obtained 
in analogy with (17) and (18a), (18b). From the condition of constancy of the incident flow 


Po. j 3 Ps u Y, (34) 


and, from the condition of constancy of the momentum 


P = 2P - Pay: 


0. > 

P. R (35) 
ins (V. À u l, eps 
" 
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8, The open boundary. Conditions on the boundary i =I are obtained in a manner 
analogous to the way we obtained (19a, b). In the present case, they give the following 
expressions for P and u in the fictitious cells (I+1, j) 


Prete j E Mia xis: , (36) 


for j = J, similar to (20), the boundary condition of equation (29) will be 


P -I 

S q 2 hd 

Vi, [sl > eni (- X P. (37) 
iş 581 i 


9. Boundary of the body. It is obvious from conditions (21) - (23) that, in the 
present case, the boundary conditions for equations (3) and (29) on i — lr will be 


^ ^ 
- * 


Pept hei Pp sui” Tipi! (38) 


on j = İyə the boundary conditions of equation (29) take the form 


P.. 
A ^ by / y 
U. - =-pP.. 1. 21 (39) 
9 T »j' I (2] ? 
671 ‘sft t ( əl 
where Pi : is the boundary condition of equation (4), the choice of which shall be made 


A 
later. 
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10. The axis of symmetry. In analogy with (24) the boundary condition of equation 


(29) on the boundary j = 0 will be 


^ 
v = AO 9 


Ú A 
5,0 ~ ” by) 
: 5 Pii 


where P; is the as yet undetermined boundary value of P of equation (4), 


, 0 


Pio (40) 


Thus, we have looked at two distinct systems of finite-difference equations approxi- 


mating the equations of hydrodynamics, and we have found the possible variants of the 
boundary conditions for these systems, 


It should be noted that, in the approximation of the equations of hydrodynamics by 
means of the system (3) - (5) with the appropriate boundary conditions (14) - (24), the 
momentum of the fluid along the z-axis is conserved, The flows of energy in the region 
of an empty cell and on the boundaries of the region in question have a completely de- 
fined meaning. When one uses the approximating system (3), (4), (29), the total energy 
of the fluid is also conserved, 


We can remove any remaining ambiguity in the choice of the proper variant of 
boundary conditions (from the initially equally-likely possible variants found above) if 


we also investigate the physical phenomena in the vicinity of an empty cell and at the 
boundary of the region, 


5, THE COMPUTATIONAL SCHEME 


1. Choice of initial conditions of the problem. 


Let us look at the initial conditions of the problem, They can be specified in different 


ways, depending on the process that we need to investigate. In the present article, we 
are primarily interested in investigating the method itself rather than any particular 


physical problem, Therefore, with an eye to saving machine time, we pose the problem 


in the following way: A cylinder that is motionless with respect to the medium under a 


specified pressure ə and density Pon achieves instantaneously a velocity Mes along the 


axis of symmetry. In this problem, we observe the onset, development, and stabilization 


of a shock wave, 


Figure 3. The Initial Distribution of Particles in 
the Cell: a)n=2;b)n=3,. 
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First of all, we choose the dimensions of the region of solution Z and R and the 
position of the cylinder in it, Keeping in mind the required accuracy of the solution and 
the possibilities of the computer, we then choose the Eulerian and Lagrangian networks, 
The dimensions of the Eulerian network are (I + 2) x (J + 2) with consideration of the 
fictitious cells. The particles of the Lagrangian network at time = 0 are distributed 
uniformly over the region, Therefore, it is convenient to put initially n? particles into 
each cell (i, j), wheren=1, 2, 3, ... . Examples of such an initial distribution of the 
particles are shown in Fig. 3. 


The initial conditions of the problem are as follows: a description of the particles 
(the kth particle has a definite mass m, and coordinates z(k), r(k) and a description of 


the Eulerian network [each cell (i, j) is characterized at time = 0 by the quantities 
M; ; I j ul j and Vi .]. Since the density of the incident flow p a is uniform, the 
m, will depend, in the axisymmetric case, on r9. 


A 
my, = "727. nArA zr" ; (41a) 
pk) _ 2In əs ar Where b - 1, ..., İn: (41b) 
(90 - ! az, where 1 - 1, ..., İn, (41c) 


where Az “ Z/L Ar =R/J. 


m the description of the Eulerian netvork, M. j is the sum of the masses of the 


particles in the cell (i, j). Therefore, at time - 0 
Mo -20 nazar" -+). (42) 


The specific internal energy corresponding to the cell (i, j) at time = 0 is determined 
from the equation of state of an ideal polytropic gas 


P. (43) 


PNE ERO M 
"pg. (Y-U0 
The velocities in the cells conform to an undisturbed incident flow: 


"jj = yo Vj = 0. (44) 
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Formulas (41) - (44) completely determine the initial conditions of the problem. 


2. Description of the Computational Cycle. 


1. Determination of the pressure. In the first step of each cycle, the pressure in 
each cell is determined from the initial values M; T and L j which depend on the selected 


equations of state, In the present article, the equation of state was taken in the form BT 


P- sly - 1) 


from which we get the pressure in cell (i, j): 


M: 
s, tds of m ES 45 
Pi; " Vij lj Cy D, 
where Vi j is the volume of the cell: 
7 2/. 1 
Vo = 2nAz(ar) ( -1. (46) 


We note that the pressure in an empty cell, which is found from Eq. (45) automatically 
vanishes, From these values Pi .is determined in a neighborhood of an empty cell, and 


also in the fictitious boundary celis: from one of the variants of the boundary conditions 
(14) - (24) or (33) - (40), depending on the finite-difference equation that is chosen to 
approximate the energy equation. 


2. Determination of the tentative new values ofu Uu, v, andi. From the values P. 


1, 


found in the preceding step, we determine the tentative new values ofu, ij and Yi 
, + 
in accordance with (3) and (4): 


: xaz(ar) |i - i 

“2 Mo, (Pian; 7” yər, (47) 
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It should be noted that the pressure in the fictitious cell Gp + 1, j) assumes different 


values in the calculation of u, : İr and V; dud In the first case, it is treated 
İr 


ə + 1° 
pa 
as the vertical boundary of the body and in the second as the horizontal boundary, For 
an empty cell, "i j and V; = have no meaning and can be arbitrary. For the sake of 


definiteness, they were taken equal to zero, 


From the values of uy . and V; N that we have found, we determine the velocities in 


the fictitious boundary cells in a mə with one of the variants of the boundary con- 
ditions (16) - (24) or (33) - (40), after which we determine the new tentative values of 
5 . from (5) or (29). These yield, respectively 


P.. fu. e as v. 
e- ^ 

= . ə 7 a fof oleh tul 
Lp, = IN BtnAz(ar) y 2 | 5 


Vy Y by] Az 
(49a) 
or 
ə P.. fu ue, 
İ, 1, , -8tzadan (ız int 7 Mi, 
: "İM, öz (49b) 
m an empty cell, I, . is taken equal to zero. For cells inside the body, we cannot make 


calculations from (47) - (49) because of the instability of equations (2). The tentative new 


values of E Lj Kz. n Kr, j are determined from (6) on the basis of the values found 
for u. 6 V. , and L. 
i, j 1, ) 1, J 


3. The motion of the particles, At this step, we integrate the equations of motion 
of the ‘particles. As calculations on a computer have shown, the integration does not 
require a high degree of accuracy and for it Euler's method is quite sufficient: 


k k) n” 
zo idi 


(50) 
(0) 6) " 
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where Ar = et/n' and Ar is the integration interval, The velocities Ya and Ya are 
determined by means of the velocity field u j and vi i The simplest procedure is to 
determine u_ and va only from the values of the velocities in the cell in which a particle 


is at a given instant, that is, 


gx 4 (51) 


Here, 


A 


m ( (4) 
Ar] - Dr, 


y <azi. 


) : " 
gar-1l, AZÜ— 1)-z 


If a particle is in an empty cell, it continues to move with the velocity it had in the 
neighboring cell. If the particle impinges on a rigid wall, it is reflected from it and the 
collision is assumed to be completely elastic. As a consequence of the collision, the /39 
component Y and Va corresponding to it changes sign [the corresponding component 
K (4) in (7) also changes sign ]. The axis of symmetry is now treated as a rigid wall. 
If the particle leaves the region, the integration is terminated and the characteristics 
of the particle m, z, and r are not stored in the memory of the computer. 


More accurate values ofu_ and va can be obtained by linear interpolation of the 


values of u and v in the four cells closest to the particle. In this case, if 


.] (2) ok NS ( k) 251 


-. 


then, as one can see from Fig. 4, 


(4) | 
iS aay dioi 
PIE e(0:7--u ( 3) ; (52a) 


140 


/” 
Dor rmn) -(i - ? (52b) 
7 Az 2/1” 


where 
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where 


Figure 4, The Linear Interpolation of Velocities 
of a Particle (k). 


If one or more of these four cells are empty, then the values u and v in them are 
taken equal to the corresponding values of u and v in the cell in which the particle is 
located at a given instant, or in which it was located at time = 0 if it is in an empty cell 
at the given instant, On the boundaries of the region of solution, we chose the following 


conditions: 


Ifi=0, then 
(53a) 
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Fori=1, 


ui) =U), ¡5 TAE = ÜL)? 
. e ı “7 (53b) 
Mika, for "Pl jar? 0575024722 
At the upper open boundary j = J, we have analogously 
Ja =p ejer = Vip 
" m ə - (530) 
Uu a “eb UPISLj* 


If the particle is close to the boundary of the body or close to the axis of symmetry, 
then the values of u and v in the fictitious boundary cells were chosen from the condition 
of impenetrability of the boundary, that is, from the condition that the normal component 
of the velocity vanishes on these boundaries, This condition was satisfied in the following 


way: 
Fori “inə 
i. o EET U. ür. 
gon = pag? "psi ” "pi" (53d) 
MIGSPES Mira jet? Vie tie. "iy, iri? 
for j= UL 
/41 
Bi jp : Hi,iT +1" ij Up aC On | (538) 
Lis jo E ZEE m7 iJ 
for j = 0, we have analogously 
Ao 2 41 : Vio 2 “ə. l (53f) 
T ; 7 -5 


If the particle began its motion iu the cell (i, j) and after a time öt was in the cell 


(i', 1), then the values of m, and (7) are taken from the values of M, j Kz, ; Kr, i 
and E, i and they are added to the corresponding values in the cell(i', ff), or they are 


omitted from the subsequent solution when the particle leaves the region of the solution. 


39 


After the integration of the motion of all particles in the region in question, the 
entry of new particles through the boundary i = 0 was simulated. This is done by the 
introduction of nJ particles with parameters (41a, b, c) with 7 = 1 into the computa- 
tional schemes, Each of these particles bears a mass m, into the corresponding cell 


(1, j). Also, 


The process of displacement of the particles and the redistribution of the mass, 
energy, and momentum between the cells of the Eulerian network is now terminated. 


As a result of the cycle of calculations, the values found for M; . and also ut 


, y 
vt ; Ti j determined from (8) are the initial conditions for the next cycle. For empty 
cells, these quantities were taken equal to zero. 


The length of the cycle 0t was taken to be 


| ne. (54) 


As was pointed out in [1], in calculations according to the particle-in-cell method 
does not require satisfaction of Courant's condition 


A 


max 


66 << AX, 


Furthermore, for very small ót, it is impossible to obtain a stable solution. The value 
of öt should now be chosen from the condition 


| V Pt "Ax. 


Tax 
The calculations by means of the cycles described above were made prior to the 


stabilization of the shock, As was shown in [4], this occurs in the present problem after 
a time 


Uta 7 — € , (55a) 


faz 


where d is the diameter of the cylinder, as is the velocity of sound in the incident flow, 


and Mo is the Mach number of the incident flow, The number of the cycles until the 


stabilization is 


S (55b) 


zb. 2n At VM. i 
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However, as calculations have shown, atf‘2r the stabilization of the flow, the values of 
its parameters in the cells continued to oscillate as a result of the fluctuation in the 
number of particles in the cells, Therefore, the values obtained for the individual cycle 
must not be regarded as the final result describing the stabilization flow. To smooth the 
oscillations in the values of the parameters of the flow, a simple averaging with respect 
to time was performed: 


5 > 
E üt. 

(s) $ = Sstab gl 56 

Bes = 

E Ss- S$ stab ” 


E M b T 7 are defined analogously. In the calculation of these quantities, 
, 9 ? 
the possibility arises of making an additional improvement in the convergence of the 
solution. In certain trial computations, after t tab is attained, the average values of 


u me V 7: and 1 7 were taken as the initial values for cycle s + 1, 


and v 


Let us now look at the distribution of the machine memory during the time of the 
calculations. The stored values of the parameters for each cell in the Eulerian network, 
at various stages of calculation of a single cycle, are shown in the following table: 


TABLE 
MIuv - the initial conditions of the cycle; 
MiuvP - calculation of the pressure via Eq. (45) and the specification of 
its boundary values on the vertical boundaries of the region and 
m of the body; | 
MIuvPu - calculation Us j from Eq. (47); 
MIuvPu - the specification of boundary values of the pressure on the 


horizontal boundaries of the region and of the body; 


MluvPuv - calculation of Vi . from Eq. (48); 


MluvPuv 


M'E'K 'K ! 
Z r 
M! I' ut y! 


the specification of boundary values of the velocities on all 
boundaries of the region; 

calculation of I from (5) or (29); 

calculation of E, j and K; j from (6); 


integration of the motion of the particles; 


calculation of the final result of a cycle via Eq. (8). 
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As indicated in the table, each cell of the Eulerian network is described by no more 
than seven parameters simultaneously, and each particle of the Lagrangian network is 
described by no more than three parameters. The total machine memory needed is 
approximately equal to 


(1 +-2X] + 237 +3n1] < 1I](7+3n>, 
where n is the average number of particles in a cell, 


6. RESULTS OF CALCULATIONS 


To make a comparison of the different variations in the method, we made several 
trial calculations. All the calculations were made on a BESM-6 computer. In all cases, 
we used the same region of solution (see Fig. 1) with the following parameters: R - 2m; 
Z x 1.4m; and I = 14; J = 20; i,, = 13; j, = 10. 


The parameters of the incident flow were taken with the following values: y = 1.4; 
o 70.132 kg + sec?/m5, Po = 10332 kg/m?; Vo, 7673 m/sec. These figures cor- 


respond to supersonic flow around an object with Mach = 2.03, 


The number of particles in a cell n? at time = 0 is equal to 4, They are distributed 
in the cell as indicated in Fig. 3a. Let us assume that the flow stabilizes after 
S stab ^ 100 cycles. The calculation was carried out to the 200th cycle. After the 100th 


cycle, the results were averaged over time. 


The undetermined boundary values of the pressure at the horizontal boundaries were 
taken in analogy with the vertical boundaries, that is, 


5. "5 Pi J under conditions (20) and (37); 
p. o GPL. 2 xcd under conditions (23) and (39); (57) 


Lj LL 
under conditions (24) and (40). 


In the trial solutions, the following questions were investigated: the development 
of a shock wave, the form of the stabilized shock wave in comparison with the experi- 
mental data of [9], the behavior of the field of flow after stabilization of the shock wave. 


The shock wave was determined from the break in the streamlines and from the 
discontinuity in density. From the increase in the density, one can determine the cell 
through which a shock wave passes and one can determine its position in that cell by a 
linear interpolation from the values of p in neighboring cells. Thus, if the shock wave 
is in cell (i, j), its coordinates were determined in the following way: 


r-ar[i— 2): add. ie O Pee : (58) 
* Fin a 
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Investigation of the method was done according to the following schemes: 


Scheme 1, In this computational scheme, the equation of the energy was taken in 
the form (5). For boundary conditions, we used P, u, and v, as determined by (14), 
(16a), (17), and (19) - (24). The velocity of the particle was determined from (51) 
without the linear interpolation (n' = 5 in (50)). 


The development of a shock wave is shown in Fig. 5. Calculations have shown 
that even in the first steps there arises a back flow in front of the body. Because of 
this, the wave propagates rapidly along the flow. On the boundary of the collision of 
the Io Ee and as back flows, there is a concentration of particles. The determina- 
tion of b^ and ve ) from (51) artificially retards the stabilization of the flow in front 


of the body and it leads to subsequent concentration of particles on the boundary between 
those cells in which the velocities are directed toward each other. From this boundary, 
the disturbances are distributed upward along the flow, as if bounced from a rigid wall. 
The shock wave at the 60th step leaves the region of the solution, During the same time, 
more empty cells appear in front of the body. This leads to high pressure gradients, as 
a result of which the particles achieve high velocities and, after several cycles, the 
flow pattern is completely distorted. 


Scheme 2. In this computational scheme, an attempt was made to decrease the 
concentration of particles in the region of the shock wave, keeping unchanged the other 
computational equations of scheme 1, For this, the following restrictions were imposed 
on the velocities of the particles: 


ds. Wee 20s (59) 


If, as a result of solving equations (3), and (4), one of the components of the velocity in 
the cell (i, j) turned out to be negative (for example, suppose that u. i,j < 0), it was 
taken equal to zero. Here, T, .1 was increased by an amount M, 4 a, ;/? in order to 
conserve the total energy of ilie gas in the cell. 


The results of the calculation are shown in Fig. 6, where the dashed curve indicates 
the shock wave obtained experimentally in [9]. In this case, the disturbances are propa- 


gated upward along the flow more slowly. For s = 60, the shape of the shock wave is close 


to the experimental one. Up to s — 100, the wave remains in the same region but no 
stabilization occurs. The effect of imposition of the restrictions (59) which violate the 
condition of conservation of momentum along the z-axis, begins to be felt, In this com- 
putational scheme, neither momentum nor energy is conserved, Beginning with s = 60, 
the wave is deformed. It is shifted toward the body at the axis of symmetry. The entire 
front of the shock wave is then distorted. Its position for s = 100 is shown in Fig. 6. 
With subsequent calculations, the shock wave continued to be adjacent to the body along 
the axis of symmetry whereas, along the boundary i - J, it is displaced slowly upward 
along the flow. With time, the flow pattern becomes more and more distorted, By the 
instant s = 80, the parameters of the flow behind the shock wave reach their stable 
values but the oscillations increase with the passage of time. The numbers in Fig. 6 
indicate the deviations (in percentages) of the parameters of the field of flow from their 
average values at the corresponding points of the region of the solution. 
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Figure 5, Results of the Figure 6. The results of 
Calculation in Accordance Calculation According to 
with Scheme 1: - - - - - , Scheme 2:- -- -- , Posi- 
Position of the Shock Wave tion of the Shock Wave as 
as Determined by Experi- Determined by Experiment 
ment in [9]; » (1), [9]; , (1), (2), (3), 
(2), (3) the Position of the (4) the Position of the Shock 
Shock Wave at Instants Wave for s — 20, 40, 60, 100. 


s — 20, 40, 60. 


Scheme 3, In contrast with the two preceding schemes, the total energy of the 
system is strictly conserved in the case of this scheme. This is done by taking the 
equation of the energy in the form (29) with boundary conditions (33c), (34), (36) - 


(40). 


In this case, the shock wave develops even more slowly. It attains a stable position 
only by s =80. Furthermore, just as with scheme 2, it begins to be deformed, Its 
position for s = 100 is shown in Fig. 7. However, these deformations are insignificant 
and the shock wave remains at all times in a neighborhood of the experimental curve. 
The oscillations of the parameters of the flow behind the shock wave decrease to about 
one-half, The calculation time to the instant of stability of flow is 5 minutes. 


Scheme 4, In this computational scheme, we investigated the influence of the 
method of determining uf) and v(9 . All the equations and boundary conditions were 


the same as in scheme 1 except that we made a linear interpolation of the velocity from 
(52) instead of (51) in the calculation of the displacement of the particles. As before, 


we took n' in (50) equal to 5. 


Just as in scheme 1, there again arises a back flow in front of the body. Up to the 
instant s = 100, the shock wave jets in the neighborhood of the experimental curve but, 
in contrast with the scheme 1, its position is stable. With subsequent calculations, the 
form of the shock wave is not deformed, and all the distortions that occur at s = 100 
are subsequently smoothed. 
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The results of the calculations are shown in Fig. 8. 


The position of the shock 


wave at s = 200 was determined in accordance with (58) from the values of p averaged 


over the time, The values of the parameters of the flow behind the shock wave con- 


tinue to oscillate around some average value even after stabilization. The oscillations 


of the averaged values are damped. Here, the average relative error far away from 


the body amounts to 3%, which is only half as large as in scheme 3, Directly in front 


of the body, these errors are significant, just as before (their values are shown in 


Fig. 8). 


Scheme 5. In contrast with scheme 4, the equation of energy was taken here in the 


form (29) with boundary conditions (33c), (34), and (36) - (40). 


The results of the calculations are shown in Fig. 9. With such an approximation of 


Figure 7, Scheme 3, 


- — --, Position of the 
Shock Wave as Indicated 

by Experiment in [9]; : 
(1) Position of the Shock 
Wave for s = 100; i 

(2) the Shock Wave from the 
Density Values Averaged 
with Respect to Time from 
s = 100 to s = 200. 


Figure 8, Scheme 4, 


~---, the Shock Wave as 
Calculated from Experiment 
of [9]; , (1) the Shock 
Wave for s = 100; , (2) 
the Shock Wave from the 
Averaged Values of p to 
the Instant s = 200. 


the energy equation, the shock wave sets in more slowly than in the scheme 4. By the 


instant s — 100, it still does not attain its stable position, The form of a shock wave 
found from the averaged values of the density s = 200 coincides well with the experi- 
mental one (see Fig. 9). The oscillations of the parameters of the flow between the 
shock wave and the body decrease vonsiderably. 


An attempt was made to increase the accuracy of the calculation by decreasing the 
interval of integration of the equations of motion of the particles, Figure 10 shows the 
results of the calculations according to the same scheme but with n' = 10, The relative 
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deviations of the parameters of the flow from the corresponding average values between 
the shock wave and the body decreased somewhat, As before, the position of the shock 
wave coincides well with the experimental position, The time of calculation up to the 
instant of stabilization of the flow with the nf = 5 amounted to 15 min and, for nf = 10, 


it amounted to 25 min, 
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Figure 9, The Results of Figure 10. Scheme 5, with 
Calculation According to n'-10: -----, the Shock 
Scheme 5 with n! = 5: Wave in the Experiment of 
----- , the Shock Wave in [9]; , the Shock Wave 
[9]; , (1) the Posi- from the Averaged Values of 
tion of the Shock Wave p for s = 200, 
for s = 100; , (2) 


the Position of the Shock 
Wave from the Averaged 
Values of p for s = 200, 


Scheme 6, Several calculations were made according to scheme 5 but with the 
addition of the restriction (59). 


The results of one such calculation are shown in Fig. 11. As it turned out, the 
restriction (59) was superfluous in schemes which included integration of the equations 
of motion of particles with linear interpolation of the velocity (52). Just as in schemes 
2 and 3, there is a deformation of the shock wave. Again, it is displaced toward the 
body at the axis of symmetry and it moves upward along the flow by the boundary j = J. 
The oscillations of the parameters of the flow behind the shock wave are somewhat 


stronger than in scheme 5, 


Scheme 7. Calculations in accordance with this scheme were made in accordance 
with the same equations and boundary conditions as in scheme 5. The initial number of 
particles in each cell was taken equal to n? = 9 and they were distributed in the cell as 
indicated in Fig. 3b. Here, the time ót of a single cycle was accordingly decreased and 
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the number of them increased. The flow stabilizes after S ab^ 150, and the calculation 
was carried out to Ssnal" 240, The time for calculation, in comparison with the cal- 


culation in accordance with scheme 5, increased by more than 300 %, 


Figure 11, The Results of the Figure 12, Calculation in 
Calculation in Accordance Accordance with Scheme 7: 
with Scheme 6: ----- , the 0  ----- , Experimental Posi- 
Experimental Position of the tion of the Shock Wave from 
Shock Wave from [9], — .. — [9]; — . . — the Shock Wave 
the Position of the Shock Wave for s = 150; the Shock 
for s = 100; the Shock Wave from the Values of p 
Wave from the Averaged Val- Averaged over the Time for 
ues of p for s = 200. s — 240. 


The results of the calculation are shown in Fig. 12. In comparison with the results 
of calculation via scheme 5, the improvements are insignificant. The position of the 
shock wave is again in good agreement with the experimental position. The deviations 
of the parameters of the flow decreased somewhat from their average values in the 
region between the shock wave and the open boundary i = I, and between the wave and 
the body they remained approximately the same. The time of calculation to stable flow /50 
amounted to 70 min. 


Scheme 8. Calculations according to this scheme were carried out in the following 
way. Up to the instant of stabilization of the flow, the calculations were made according 
to scheme 5. Then, for every s» s stab’ the average values of uj S), V {s 7 and Tİ >, 


1, 


in (56) taken as initial conditions foi the (s + 1) cycle. The ecuador and Pounders con- 
ditions remained analogous to scheme 5. 


The results of the calculations are shown in Fig. 13. The position of the shock wave 
is in good agreement with the experimental position. After s = 100, the solution rapidly 
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converges. After 40 cycles, after stabilization of the 
flow, the deviation of the parameters from their 
averaged values does not exceed 1% over the entire 
region behind the shock wave except in a neighborhood 
of the stagnation point. Here, it should be noted that 
a significant increase in the accuracy for calculations 
in accordance with scheme 8 is achieved with hardly 
any increase in the calculating time, Figure 13 shows 
the streamlines corresponding to the time s = 200. 
The time of calculation to the instant of stable flow 
amounted to 10 min, and the total time of calculation 
amounted to 25 min. 


To illustrate this scheme, we calculated the flow 
around a cylinder for Mo = 2.91 and Mo = 4,35, In 


the first case, we set Y = 966.3 m/sec, ip = 10, 


Figure 13, Scheme 8: 


Mə - 2.08: ----- , the B db 200, Sánal ” 320; in the second case, 
Shock Wave from Experi- V_ = 1444,0 m/sec and i, = 9, 
oo T 

ment [9]; —. . —, the 
Shock 5: m 7 7 The remaining initial and boundary conditions 
2 7 ə are analogous to the preceding case, The positions 

d Values of o for of the shock wave and the acoustic line obtained in 
7 2 th the calculation according to scheme 8 for the given 
EE 0 values of the Mach numhber are shovn in Fig. 14, 


Streamlines Correspond- 
ing to the Instant s — 200, 


2,0 LEE ao E 


Figure 14. a) Mə = 2.91, b) Mo —4,35, The Results of the Cal- 


culation According to Scheme 8: ----- the Shock Wave from the 
Experiment of [9]; the Calculated Shock Wave for s = 320; 
—.. — the Acoustic Line. 
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7. CONCLUSION 


The calculations according to all the schemes yield Py i =P eae The boundary 


conditions (18a, b), (17) and (34), (35) coincide; that is, for supersonic flow, the con- 
ditions on the boundary i = 0 have the form (17) or (34). 


It is desirable that the value of the pressure be specified on the open boundaries 
more accurately since conditions (19a) and (36) yield values whichare too high, and this dis- 
torts the flow in the neighborhood of the open boundaries. Therefore, the acoustic line 
left the region of the solution in all the calculations. 


In comparing the computational schemes, we can draw the following conclusions: 


1. Schemes 1 and 2 are not applicable for analysis of flows with great deformations. 


2. Scheme 3 can be applied in the case of tentative calculations of nonstationary 
flows if great accuracy is not required. In the stationary case, it is sensible to carry 
out the solution only up to the instant of stabilization of the flow. This scheme can be 
used for nonstationary flows without great deformations when a slow computer with a 
sufficiently large memory is employed. 


3. With calculations of nonstationary flow that require a specified degree of 
accuracy, we can use scheme 5, The decrease in the integration interval of the equa- 
tions of motion of the particles produces an insignificant increase in the accuracy but 
a large increase in machine time. In this scheme, it is sufficient to choose n' between 
3 - 6. In calculations of the flows involving subsequent compression of the medium, 
increase in the initial number of particles in a cell does not lead to a great increase in 
accuracy although, as one can see from the calculations, it does considerably increase 
machine time. The accuracy of calculations according to this scheme can be increased 
noticeably by decreasing the dimensions of the cell of the Eulerian network. In the case 
of calculations of flows with a fine structure or of region of special interest (for ex- 
ample, the neighborhood of a body), we can use locally a finer Eulerian network. Here, 
it is necessary to partition a particle entering this region into several smaller ones. 

It is probable that this will lead to an increase in the computational accuracy. 


4, For calculation of stationary axisymmetric flows, we can use scheme 8. All 


the remarks about'the use of scheme 5 hold for this scheme. It should be noted that, for 


a more rapid convergence of the solution, it is necessary to determine the instant of 
stabilization of the flow as accurately as possible, 


REFERENCES 


1. Harlow, F. H. Hydrodynamic problems involving large fluid distortions, Journal 
of the Association for Computing Machinery, April 1957, Vol. 4, No. 2. 

2, Buneman, O. Dissipation of Currents in Tonized Media, Phys. Rev., 1959, 
115, No. 3. 


3. Lomnev, S, P. Influence of Coulomb repulsion on the grouping of particles in a 
linear electron accelerator, Doklady AN SSSR, 1960, 135, No. 4. 
4, Evans, M. W., and F. H. Harlow. Calculation of supersonic flow past an axially 


symmetric cylinder, Journal of the Aeronautical Sciences, April 1958, Vol. 25, 
No. 4. 


/52 


49 


© 00 


Evans, M. W., and F. H. Harlow. Calculation of Unsteady Supersonic Flow Past a 
Circular Cylinder, ARS Journal, January 1959, Vol. 29, No. 1. 

Evans, M. W., F.H. Harlow and B. D, Meixner. Interaction of Shock or Rarefaction 
with a Bubble, the Physics of Fluids, June 1962, Vol. 5, No. 6. 

Harlow, F.H, The Particle-in-Cell Method for Numerical Solution of Problems 
in Fluid Dynamics, Experimental Arithmetics, High-Speed Computing and Mathe- 
matics, 1963, Vol. 15. 

Computational Methods in Hydrodynamics, Moscow, Mir Press, 1967. 

Maslennikov, V.G. Onthe Form of a Leaving Shock Wave Formed in the Case of 
Supersonic Motion of a Hemisphere and a Cylindrical Block in Various Gases; In 
the Collection: Aerofizicheskiye issledovaniya sverkhzvukovykh techeniy (Aero- 
physical Investigations of Supersonic Flows), Moscow, Nauka Press, 1967. 


/53 


SOME QUESTIONS IN THE CALCULATION 
OF THE EVOLUTION OF STARS 


S. N. Lomnev, K. A. Begishkina, 
Z. F. Levina, G. V. Ruben 


ABSTRACT: Application of the hydrodynamic PIC method to 
the problem of evolution of stars is discussed. 


For the past several years there has appeared in the literature a large number of 
articles devoted to calculation of the internal structure and gradual development of 
one-dimensional models of stars. The problem is formulated as a Cauchy problem for 
the system of equations 


ax, t): 
3i f (P, R, L, T, 0 (1) 


oY . P 5 t , 
at PA ,R,L T, ) 


describing the processes of nuclear transmutations of hydrogen and helium over time, 


The internal structure of a star for fixed t is determined from the solution of the 
boundary value problem for the system of nonlinear ordinary differential equations 


de ^ 

dP -4 

27: (2) 
qe CLIC, X* pr* XX VPO T ^1; 

dq dinP, dq i 


where the Ci are constants and p is the density inside the star, which relate the varia- 


tion of the radius R, the pressure P, the temperature T and the luminosity L with 
respect to the mass q of the star. 
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The basic difficulty is embodied in the solution of the boundary value problem, /99 
Earlier presented articles have been based on two methods: "trial integrations" 
(described by Schwartzschild [1]) and the double-sweep (proposed by Henyey, Levier 


et al, [2]). 


The first uses a well known approach: the combination of trial integrals of (2) from 
two sides (from the surface and from the center of the model) with approximately given 
"missing' boundary conditions and calculation of corrections to these approximately 
given boundary conditions, To the disadvantages of the method one must ascribe the 
poor convergence of the corrections, Under the burning out of hydrogen in the center 


it is useless. 


The double-sweel method offers the possibility of calculating the evolution of stars 
up to the late stages of burning of hydrogen and the onset of the burning of helium. How- 
ever, using the linear approximation of the equations, we obtain a more precise value 
of the calculated quantities. 


The difficulties of solution of two-dimensional problems and nonstationary processes 
of stars by these methods oblige one to seek new paths of research. 


One of the possible processes of calculation of two-dimensional evolutional problems 
can be the following. The system (2) is split into two subsystems, which are solved suc- 
cessively. In the capacity of the two "missing" functions are taken the values of the pre- 


ceding model. 


The first pair of equations will be expressions for the variation of P and R, the 
second pair for T and L, Such a choice of pairs of functions is illuminated by the fact 
that P or and R characterize matter, and T, L characterize the fields created by this 
matter. Thus, matter evokes variation of the fields, and the fields variation of matter. 
For a sufficiently fine step relative to t one should expect a sufficiently precise solution. 


The splitting of (2) offers the possibility of using several schemes of calculation 
relative to t, The following appears interesting to us: 


Scheme 1. Ateach step At one solves first the system for P and R. Then, not 
changing t, with the newly calculated R4 (a) and P4 (q) the system for L and T is inte- 


grated, The newly calculated L4 (q) and Ty (q) serve for the determination of R, (q) 
and Po (q) up until the (i + 1)-th expressions do not differ from the i-th within a pre- 


scribed small magnitude. 


Scheme 2, Pairs of equations are evenly divided with respect to time: in the odd 
steps R, P are calculated, in the even L, T. 


Scheme 3. For fixed t the system for R and P is first solved. The calculated 
values are used in the integration of the equations for L and T. 


Test calculations of the evolution of one-dimensional models stars were carried 29 
out. The best results were obtained with the use of scheme 2, Thus, on the time interval 
0-10 years agreement was obtained with the magnitudes, calculated by another method, 


of 4-5 symbols, 
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Utilization of splitting of the equations appears advantageous to us in the later 
stages of evolution, for example, in time intervals where the Schwartzschild method 


converges poorly. 


One-dimensional stationary processes with scales of phenomena of 109 give certain 
averaged characteristics of the development of stars, Hydrodynamic instabilities and 
two-dimensional phenomena (with scales of about a year and less) are more advanta- 
geously interpreted by means of ''coarse particles’. This offers a series of advantages. 


1. The system of equations is simplified and is decomposed into two groups of 
equations which are sequentially solvable, 


2. The equation of continuity, which is unstable for numerical solution, drops out 
and is replaced by a Cauchy problem for a system of ordinary differential equations. 


Strong deformations of matter are permitted and, naturally, the physical process 
is described, 


In this paper two versions of the method of "coarse particles" are employed: the 
method of "mean density" [3] and the method of "particles-in-cells" [4] (PIC for short). 
For exhibiting the particulars of the calculation the method is tested on a spherically 


symmetric star. 


The "mean density" method, The radius of the star is divided into N invervals 
(cells with indices i - 1, 2, ..., N), The mass of each cell is decomposed into / 
particles. The calculation is begun by the determination of the density drom the dis- 
tribution among the cells. From the equation of state for the density the pressure 


C oT 
P. = 15 


m Ə 
ç 9 + +(C,1)', (3) 


is obtained, where yu is the mass of the i-th molecule, C1; Cy are constants; the 


quantities p i and Pi relate to the center of the i-th cell, 


The next step consists in the integration of the field equations 


- 2xr7 ple - Er); 
o (4) 


qo 
ww 


aT | dinT , GMe 
ar | din P r” i 


where £ is the energy of burning, 8 G is the energy of contraction or expansion, 


r 
M = .. 4T r?pdV, Gis the gravitational constant, with the boundary conditions 
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L(0,t) = 0 


, ~ (5) 
T(r, ,t = L (rr Dr. T. , 


where T eff. is a constant, 


The system (4), (5) is solved with the aid of trial integrations from both ends and 
determination of corrections at the point of joining. In the calculated fields motion of 
the particles takes place according to the equations 


dt ” dt P or TM (6) 


The system (6) was integrated by Euler's method. The gravitational force GM/r? 


was determined relative to particles, the pressure gradient I p relative to cells. 
The form of the initial values us (0) and r. ¿(0) depends on the “. problem. 


"Thus, calculation of the evolution of a star on the time interval T= .. At ,, con- 
m 


sists of calculation of the field from (4), (5) for fixed t and shifting the masses (6) in 
these fields by At wherein At is selected in accord with the condition of prohibition 


of transition of particles between cells ór, = u, At 


The distribution of particles and cells appears to be a complex problem and depends 
on a series of causes, Substantial difficulty resides in the sharp drop in density and 
pressure, changing by approximately 10 orders of magnitude from the center to the 
surface, 


In this article two forms of the distributions of particles and cells were used: 


a) the particles were of the same mass and were distributed proportionally to the 
initial density, The cells were taken such that into the first one four particles fell, 
Thus, on the surface and at the center of the star large cells were obtained; 


b) particles and cells were distributed from the condition of small relative varia- 
tions of the quantities of the field and pressure of the initial uniform model; 2000 par- 
ticles of various masses were placed in 200 cells (10 in each). Thus, cells were 
obtained which were small at the center, increased toward the middle portion, and 
then again decreased toward the surface, 
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These two forms of the distributions were used in the calculations, The second /60 
type of distribution of particles and cells gave more stable results of the calculations, 


A series of general deficiencies were revealed. 


1. In view of the constancy of the pressure gradient within a cell there are no 
collisions between particles, Hence particles of a cell move together and coagulation 
is obtained. (Fig. 1). 


2. Smalldisplacements of particles produce large variations of their velocity, 
and at each step they are accumulated. Attempts to smooth the pressure gradient and 
thereby decrease the scatter of the quantities did not permit complete removal of com- 
putational instabilities. 


The method of "particles in cells! [4]. "The deficiencies of the "mean density" 
method necessitated turning to the method of "particles in cells", which is used in 
certain complex nonstationary problems of hydrodynamics. Matter is divided into a 
large number of "coarse particles" (Lagrangean mesh), which move like points, having 
a given mass, across cells (Eulerian mesh), The density of matter in a cellis equal to 
the ratio of the mass of all the particles found in it to its volume, It is constant within 
a cell and changes stepwise between neighboring cells, The field quantities (i.e. , tempera- 
ture, energy flow, etc.) are also constant within a cell and are determined from difference 
equations, Solution of the problem leads to a series of steps according to time, Each step 
consists of three stages: calculation of the field quantities, movement of the particles 
and, finally, corrections of the field. 


Calculations carried out by Harlow and others have shown that the method is ex- 
pedient for problems in which large fluid deformations and large relative displacements 
are encountered, and impingement of the surfaces of separation of media arises. 


In the problem considered there must not be regions of flow with small subsonic 
velocities. Inside large regions of flow there must not be encountered small regions in 
which it is necessary to know the solution in detail. We undertake to apply the scheme 
of the PIC method for the calculation of fast processes in the interiors of stars. 


We shall consider the adiabatic case. We separate the star into N spherical zones 
identical in radial dimension (Fig. 2). In the first stage the fields are calculated accord- 


ing to the equations 
u.—u. M.m. 
i Ü d 2 
i At J p2 ( I. m 


i 


(7) 
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here m, is the mass of the i-th cell, Uu. is the velocity, Pi the pressure, Qi the internal 
energy, and Vi = 4mri ÖT, is the volume of the i-th cell, and “ denotes the results 


after the first stage 
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The boundary conditions are obtained from the laws of conservation of energy 


i 0; P = 0; 


N+1 7 N +1 
22[P, (1), , «P, , 01700, 1=0, 


where N is the last cell on the surface. 


Forn=1, i.e. atthe center, this equality is satisfied out of symmetry. In empty 
cells the pressure and velocity are equal to zero, 


Initial conditions are taken from the results of calculation of the initial uniform 
model, The star was divided into 200 cells identical in dimension which were fixed, and 
to them were added another 100 empty cells, In each cell were distributed uniformly 10 
particles of like mass for a given cell. 


Partition into identical cells is advantageous for calculation and insures fulfillment 
of the law of conservation of energy, although it is crude for the surface, and in a 
central cell in the calculation of volume by the formula Vi = 4T Ti Or, the relative error 


amounts to 25%. 


In the second stage particles move with velocities obtained by linear interpolation 
on the centers of two neighboring cells 


here Ty ry are the coordinates of the particle before and after the movement, Vk is 


the velocity of the particle. 


Trial calculations showed significant improvement of stability of the scheme in 
comparison with the case where the particle velocity was taken equal to the value at 
the center of the cell, 


The Eulerian quantities are not altered if a particle does not transgress the 
boundary of a cell, In the contrary case the change is determined by addition and 
subtraction of masses, energies (kinetic heat advantage) and the amount of movement 
with the corresponding quantities by a given particle, The amount of movement of the 
particles is Y, = m,U,, the energy is E, = Em /n,, where mə U, and E, are the 


parameters of the old cell after the first stage. In the third stage the corrections of the 
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field are calculated. After the movement of the particles we have Yj, Ej, mj as the 


impulse, energy and mass of the i-th cell, We calculate 


ə e .2 


Y. , . M. 
n=; Q.-E.- ———; T,, P,, e 
wm 


” 
$ 


for the i-th cell, These values are taken as the initial data for the next step of the 
calculations, 


The masses of particles at the center and on the surface differ by 7 orders of 
magnitude. The feasibilities of the machine do not permit smoothing this difference 
completely; hence the calculations have a trial methodological character. An attempt 
to calculate pulsations of stars did not give good results, the method is too coarse for 
small pulsation velocities. The errors arising under sufficiently large time of compu- 
tation exceed by far the given initial velocities. In the equilibrium model under con- 
sideration with zero initial velocities a velocity of movement arises and increases as 
a result of errors. In central cells the parameters are agitated in the first steps of 
the calculation. It is necessary to decrease substantially the time step or to throw out 
the n central cells. Discarding the center, the boundary conditions 


were introduced. Particles are elastically reflected from the boundary of cells n and 
n+1, With the increase of n the time step required for the calculation increases, for 
n=50, At =0,0005 hours, for n = 155, At = 0.005 hours, For n = 50 after 153 steps 
it turns out that the greatest velocities occur in the central cells, at the surface, and 
in the neighborhood of ri» Yo (Fig. 3). Up to this moment of the calculation the density 


has not changed, except the density of certain cells at the center and at the surface. 
Thus, the velocity characterizes instability in the first step of the calculation. 


An attempt to improve the accuracy by breaking down the cells for that same 
quantity of particles (they took N = 400 and 2000 particles) led to poorer results 
(Fig. 4). The pattern is particularly unstable in the neighborhood of Ro: In the re- 


maining cells the results do not differ significantly in absolute value from the results 
for N = 200. In Fig. 5 are presented the values of the velocity for T =0.07 hours and 
for T =0.27 hours. The calculations were carried out with a step of 0.005 hours, half 
of the radius of the star being fixed. It is evident that the most stable part is the vicinity 
of Ry. Forr > 3R. p/4 the increase of velocity is less significant. 


59 


/63_ 


AnH 


Particles 


——_ hz Rogp/00. 


Let us introduce the following distribution of velocity: u = 6R y /T for /6 
osrés Ro « u=0 for r > Ro: This is equivalent to the formation of a shock wave. 


Calculations were performed with Ro = RosB/8 (T — 1 hour: Ro is the radius of 
the sun ): t = 0.0001 hours. 


The initial P, p, T were taken from the results of the calculation of the uniform 
model, The calculations showed that the shock wave front diffuses with the course of 
time (Fig. 8), which indicates the existence of viscosity, not taken into account by the 
usual PIC method. 


One of the defects of the method consists in the crossover of particles. This 
phenomenon can be averted in the following way. Let v, r be the velocity and the 
coordinate of a particle; v is interpolated relative to the cells 


UR. ya: (8) 
3 r. - Tf, : 
bol i 
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In the interval At the velocities of the cells do not change, At time t , r-r.. 


Solution of (8) gives 


ita 6 (u T.—4.T. ,) isi 34 


(t)-e pr PA 2 Bcc d 
u. —— u . 
ixl ¿ 


If a particle crosses the middle of the i-th cell at time t' = İm + Af", and falls into 


interval 2 (see Fig. 2), then 


Particles in the inverval At do not cross over. One must turn attention to artificial 
viscosity in equations (7) and to a rational distribution of particles and cells (for 
example, a combination of a and b). The instabilities portrayed in Fig. 4 can be 
significantly decreased, 


In conclusion it should be noted that, despite the limited possibilities of the method 
of "particles in celis" for the calculation of the evolutions of stars, as it appears to us, 
this methodology is applicable to the analysis of explosive processes. 
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